TD1 - Econometrics

Partie 1 : Prise en main des dates, ts object, charts, loops and functions

a - Time series class object

Génération de series temporelles en R :
  • Une série temporelle peut être perçue comme une liste de nombres comprenant les temps auxquels ces nombres ont étés enregistrés. Dans cette première partie, on s’intéresse à l’information qui peut être stockée dans les objets de type “ts” (Times series) sur R.
    • En pratique : un objet ts(vector, start=, end=, frequency=) convertit un “vecteur numérique” en objet “série numérique”" R.
  • start/end : temps de la première et de la dernière observation.
  • frequency : le nombre d’observations par unités de temps (1=annual, 4=quartly, 12=monthly, etc.). 
    • tseries <- ts( rnorm (100) , start = c(2000 , 1) , frequency = 4)
print ( tseries )
    ##              Qtr1         Qtr2         Qtr3         Qtr4
## 2000 -0.114512864  0.371500957 -1.050950581  0.173428169
## 2001  1.209290740 -0.844565861 -0.274170303 -0.133457253
## 2002  0.110610663  2.083866977 -0.014200290  2.345218433
## 2003  0.437083992 -1.933461362  0.079837759 -0.908015208
## 2004 -0.939407840 -1.341366698 -0.151374817 -0.991097612
## 2005 -0.092510225  1.546928731 -2.273574890  1.157886377
## 2006  0.308194975 -0.570808879 -0.245915396  0.379246433
## 2007  0.550353428 -0.777396598 -0.053898443 -0.008679362
## 2008 -0.118543293  0.691106777 -0.075377029  0.355871024
## 2009 -1.165070153 -0.124902106 -1.904143703 -1.224714602
## 2010 -1.223679851  0.704920391 -0.944621779 -0.230030273
## 2011  0.962365888  0.829706516 -2.117548967 -1.033408317
## 2012 -0.217051691 -1.078160300 -0.796513101  1.156743830
## 2013 -1.507012714  0.572066880 -0.832517827  0.188429141
## 2014 -0.198565961 -0.778030929  0.496183652 -0.320770848
## 2015  0.253771132 -1.353279445  0.010929813  0.524078095
## 2016 -1.766937009 -0.251877826  1.817532838  0.621052724
## 2017  0.701873498 -0.660270934 -0.143302121  0.427095807
## 2018  0.449196392  0.598740622 -2.029406151 -0.173749103
## 2019  0.869157856 -1.255830753  0.469166919 -0.720760905
## 2020 -2.034796552 -0.276499343 -1.706079486 -1.011353714
## 2021 -0.347743654  0.493544122 -0.141839855  0.328228755
## 2022 -2.255928460 -1.300906570  0.798050555 -0.795941709
## 2023  0.180852712  0.079990577 -1.918448417  0.427797017
## 2024 -0.936514549 -0.450056454 -0.296489714  0.871835518
    #On obtient un affichage trimestriel. le vecteur Rnorm passé compte 100 valeurs. à frequence trimestrielle, celà correspond à 25 ans. Le résultat est cohérent.
    Si l’on veut réduire la donnée à l’information comprise entre le premier trimèstre de l’an 2000 et le dernier de l’an 2012, on utilise la fonction window :
  • window(x, start = NULL, end = NULL, frequency = NULL) est une fonction générique qui extrait la une sous chaine d’un objet x entre la date start et la date end, elle peut être réarrangée différemment de l’objet de départ en choisissant une nouvelle frequence. 
    • tseries_sub <- window ( tseries , start =c(2000 , 1) , end=c (2012 ,4))
print ( tseries_sub)
    ##              Qtr1         Qtr2         Qtr3         Qtr4
## 2000 -0.114512864  0.371500957 -1.050950581  0.173428169
## 2001  1.209290740 -0.844565861 -0.274170303 -0.133457253
## 2002  0.110610663  2.083866977 -0.014200290  2.345218433
## 2003  0.437083992 -1.933461362  0.079837759 -0.908015208
## 2004 -0.939407840 -1.341366698 -0.151374817 -0.991097612
## 2005 -0.092510225  1.546928731 -2.273574890  1.157886377
## 2006  0.308194975 -0.570808879 -0.245915396  0.379246433
## 2007  0.550353428 -0.777396598 -0.053898443 -0.008679362
## 2008 -0.118543293  0.691106777 -0.075377029  0.355871024
## 2009 -1.165070153 -0.124902106 -1.904143703 -1.224714602
## 2010 -1.223679851  0.704920391 -0.944621779 -0.230030273
## 2011  0.962365888  0.829706516 -2.117548967 -1.033408317
## 2012 -0.217051691 -1.078160300 -0.796513101  1.156743830

    b - Multiple time series

    On peut tracer plusieurs séries temporelles en les combinant dans une matrice.

      # Get the data points in form of a R vector .
ts.1 <- c(799 ,1174.8 ,865.1 ,1334.6 ,635.4 ,918.5 ,685.5 ,998.6 ,784.2 ,985 ,882.8 ,1071)
ts.2 <- c(655 ,1306.9 ,1323.4 ,1172.2 ,562.2 ,824 ,822.4 ,1265.5 ,799.6 ,1105.6 ,1106.7 ,1337.8)

  # Convert them to a matrix .
matrice <- matrix (c(ts.1 , ts.2) , nrow = 12)
matrice
    ##         [,1]   [,2]
##  [1,]  799.0  655.0
##  [2,] 1174.8 1306.9
##  [3,]  865.1 1323.4
##  [4,] 1334.6 1172.2
##  [5,]  635.4  562.2
##  [6,]  918.5  824.0
##  [7,]  685.5  822.4
##  [8,]  998.6 1265.5
##  [9,]  784.2  799.6
## [10,]  985.0 1105.6
## [11,]  882.8 1106.7
## [12,] 1071.0 1337.8
    # On indique le nombre de ligne que l'on veut et la répartition est faire pour respecter la dimension.

  # Convert it to a time series object .
rainfall.ts <- ts( matrice , start = c (2012 ,1) , frequency = 12)
rainfall.ts
    ##          Series 1 Series 2
## Jan 2012    799.0    655.0
## Feb 2012   1174.8   1306.9
## Mar 2012    865.1   1323.4
## Apr 2012   1334.6   1172.2
## May 2012    635.4    562.2
## Jun 2012    918.5    824.0
## Jul 2012    685.5    822.4
## Aug 2012    998.6   1265.5
## Sep 2012    784.2    799.6
## Oct 2012    985.0   1105.6
## Nov 2012    882.8   1106.7
## Dec 2012   1071.0   1337.8

    Et utilisation des fonctions d’union/d’intersection de séries temporelles

    ts1 <- ts(101:106 , freq =4 , start =2012.75)
ts1
    ##      Qtr1 Qtr2 Qtr3 Qtr4
## 2012                 101
## 2013  102  103  104  105
## 2014  106
    #intesection : affiche les valeurs interséctées, à ou il y'a des NA, rien n'est affiché
lag(ts1,-1)
    ##      Qtr1 Qtr2 Qtr3 Qtr4
## 2013  101  102  103  104
## 2014  105  106
    ts.intersect(ts1 , lag( ts1 , -1)) 
    ##         ts1 lag(ts1, -1)
## 2013 Q1 102          101
## 2013 Q2 103          102
## 2013 Q3 104          103
## 2013 Q4 105          104
## 2014 Q1 106          105
    #union : affiche tout
ts.union(ts1,lag(ts1,-1))
    ##         ts1 lag(ts1, -1)
## 2012 Q4 101           NA
## 2013 Q1 102          101
## 2013 Q2 103          102
## 2013 Q3 104          103
## 2013 Q4 105          104
## 2014 Q1 106          105
## 2014 Q2  NA          106
    cbind ( ts1 , lag ( ts1 , -1))
    ##         ts1 lag(ts1, -1)
## 2012 Q4 101           NA
## 2013 Q1 102          101
## 2013 Q2 103          102
## 2013 Q3 104          103
## 2013 Q4 105          104
## 2014 Q1 106          105
## 2014 Q2  NA          106

    c - Date class object

    Découverte des différentes options pour gérer les dates en R

    #Creation d'une date
as.Date('1915-6-16')
    ## [1] "1915-06-16"
    as.Date('1990/02/17')
    ## [1] "1990-02-17"
    #les deux formats ci dessus sont reconnus par défaut.

#On peut également choisir de préciser notre propre format : 
as.Date('1/15/2001',format='%m/%d/%Y')
    ## [1] "2001-01-15"
    as.Date('April 26 2001',format='%B %d %Y')
    ## [1] NA
    period <- as.Date (c("2007-06-22", "2004-02-13"))
  #On peut également calculer la durée entre deux dates:
days <- period[1] - period[2]
days #Durée entre les deux dates passées en paramètres.
    ## Time difference of 1225 days

    d - Plotting Time series

    Il existe deux fonctions principales pour tracer des séries temporelles. (ts_ plot : plus simple, ggplot2 : plus riche).

    Package TSstudio : starting with ts_plot

    #install.packages ("TSstudio")
library ( TSstudio )
data ( USgas )
ts_info ( USgas )
    ##  The USgas series is a ts object with 1 variable and 238 observations
##  Frequency: 12 
##  Start time: 2000 1 
##  End time: 2019 10
    #Tracé
ts_plot ( USgas ,
           title = "US Monthly Natural Gas Consumption",
           Xtitle = "Time",
           Ytitle = "Billion Cubic Feet",
           slider = TRUE)
    # slider permet d'ajouter un curseur en bas du tracé qui permettra de personnaliser la longueur de la fenêtre du tracé
    data ("Coffee_Prices")
#ts_info( Coffee_Prices )
ts_plot ( Coffee_Prices , type = "multiple")

    Package ggplot2:

    # Libraries
library ( ggplot2 )
library ( dplyr )
# Generating data within an uniform distribution
data <- data.frame (
  day = as.Date ("2017-06-14") - 0:364 ,
  value = runif (365) + seq ( -140 , 224)^2 / 10000)

# Most basic bubble plot
p <- ggplot(data , aes ( x = day , y = value )) + geom_line () +  xlab ("")
p + scale_x_date ( date_labels = "%b")

    #Autres exemples de personnalisation : 
#p + scale_x_date ( date_labels = "%Y␣%b␣%d")
#p + scale_x_date ( date_labels = "%W")
#p + scale_x_date ( date_labels = "%m -%Y")

# Libraries
library (ggplot2)
library (dplyr)
library (hrbrthemes)
# Dummy data
data <- data.frame (
  day = as.Date ("2017-06-14") - 0:364 ,
  value = runif (365) + seq ( -140 , 224)^2 / 10000)
# Most basic bubble plot
p <- ggplot (data , aes ( x = day , y = value )) +
  geom_line ( color ="steelblue") +
  geom_point () +
  xlab ("") +
  theme ( axis.text.x = element_text ( angle =60 , hjust =1)) +
  scale_x_date ( limit =c(as.Date ("2017-01-01") ,as.Date ("2017-02-11" ))) +
  ylim (0 ,1.5)
p

    #### e - Loops and conditional expressions : Boucles F and FOR sur R

    a=0
if(a!= 0){print (1/a)
}else{ 
print("No reciprocal for 0.")}

for ( i in 1: 4){
  print(log( i ))}
for ( i in c( -8 , 9 , 11 , 45)){
  print ( i^2)}

fruits <- list ("apple", "banana", "cherry")
for ( x in fruits ) {
  if ( x == "cherry") {
    break}
  print ( x )}

for ( i in 1:3){
  for ( j in 1: i ){
    print ( i + j )}}
    • Application
    mt = 5
nt =5
cter =0
mym = matrix (0 , mt , nt )
for ( i in 0: mt ){
  for ( j in 0: nt ){
    if( i == j ){
      break ;}else {
      mym [i , j ] = i*j
      cter = cter +1}}
  print ( i*j )}
print ( cter )

#### f - Next and Break Statement
x <- 1:5
for ( i in x ){
  if ( i == 3){
    break}
  print ( i )}

x <- 1:5
for ( i in x ){
  if ( i == 2){
    next}
  print ( i )}

    Partie 2 - Exercices

    La partie qui suit consistera
  • Dans un premier temps: à modéliser le comportement d’un bruit blanc gaussien ainsi que l’allure de son graphique d’autocorrélation.
  • Dans un second temps, à étudier le comportement d’une marche aléatoire (somme de bruits blancs gaussiens).
  • Enfin, à appliquer aux deux modèles le test de Ljung-box pour determiner de la manière la plus quantitative possible l’existence ou non d’autocorrélation.

    • Exercice 1 - Etude d’un bruit blanc

    1 - a - Studying a Gaussian White Noise

  • Generate 1000 observations from a normal distribution with 0 mean and a standard deviation of 1.7. This process is named (ut).
  • Using a function plot, plot the generated series.
  • Explain what is an autocorrelation function. Plot the ACF of (ut). How do you interpret the ACF results ? Are they aligned with the definition of the white noise. 
    • library(ggplot2)
ut <- rnorm(1000,0,1.7) #1000 observations d'une distribution normale de moyenne nulle et d'écart type 1.7
dataf <- data.frame(x = 1:1000,y = ut)
ggplot(dataf, aes(x = dataf$x, y = dataf$y, fill ="legend")) + geom_line(color = "steelblue") + xlab("Index") + ylab("ut") + ggtitle("White Noise ut")

    L’autocorrélation (ou l’autocovariance) d’une série fait référence au fait que dans une série temporelle ou spatiale, la mesure d’un phénomène à un instant t peut être corrélée aux mesures précédentes (au temps t−1, t−2, t−3, etc.) ou aux mesures suivantes (à t+1, t+2, t+3, …). Une série autocorrélée est ainsi corrélée à elle-même, avec un décalage (lag) donné.

  • acf : The function acf computes (and by default plots) estimates of the autocovariance or autocorrelation function. Function pacf is the function used for the partial autocorrelations. Function ccf computes the cross-correlation or cross-covariance of two univariate series.
  • plot.acf : Plot method for objects of class “acf”. 
    • Autocorrelation <- acf(ut, lag = 5)
    On remarque uniquement la présence d’un pic en 0 (qui corréspond au temps t) alors que toutes les autres autocorrélations sont significativement nulles (statistiquement) à partir de t-1. On peut ainsi dire que le procéssus n’est pas autocorrélé du tout. De plus, le resultat observé s’aligne avec la définition du bruit blanc gaussian qui est un “processus sans mémoire”.

    1 - b - Studying a random walk

    1. Generate a random walk based on the simulated white noise (keep the same seed).

    2. Convert the simulated data into a ts object (choose the the start date and freq). Plot the generated series. Plot its density. What do you observe ?

    3. Plot the ACF of (wnt). How do you interpret this results ? Can the models seen during the class be used for this type of data ? Why ?

    4. Run the Ljung Box test for the processes given in the exercise 1-a and 1-b and for a different lags (a for loop could be usefull). Plot the values of the Q(m) and conclude

    Une marche aléatoire est une somme de bruits blancs gaussiens. On a déjà définit précedemment le bruit blanc gaussien “ut”
  • cumsum() Returns a vector whose elements are the cumulative sums, products, minima or maxima of the elements of the argument.
    • En pratique, cumsum(ut) renvoie la somme de toutes les valeurs de ut (ie : ut + ut-1 … ut-k). 
    #1.
yt <- cumsum(ut)

#2.
yt.ts <- ts(yt, start = c(2000,1), frequency = 12)
plot(yt.ts, type = "l", main = "Random Walk")

    plot(density(yt.ts))

    On constate que la représentation graphique de la densité est celle d’une loi normale. Ce résultat est cohérent dans la mesure où la marche aléatoire générée est une somme de bruits blancs gaussiens. De plus, la marché aléatoire générée s’apparente à un processus à moyenne mobile de paramètres tous égaux à 1. 
    acf(yt.ts, lag = 20)

    L’acf montre que le processus est entièrement autocorrélé. Avec les processus MA et AR, on ne peut pas déduire l’ordre d’un tel acf car toutes les autocorrélations sont statistiquement différentes de 0.

    Test de Ljung-Box

    Etant donnée que l’on travaille avec une série à moyenne constante, on peut appliquer le test de Ljung-Box pour établir la présence d’autocorrélation ou non dans le processus stochastique.

  • Box.test(x, lag = 1, type = c(“Box-Pierce”, “Ljung-Box”), fitdf = 0) : Compute the Box–Pierce or Ljung–Box test statistic for examining the null hypothesis of independence in a given time series. These are sometimes known as ‘portmanteau’ tests.

    • Test de Ljung-Box appliqué aux résultats des parties 1-a et 1-b :

    ut.test <- Box.test(ut,lag = 20, type =  "Ljung")
yt.test <- Box.test(ut, lag = 20, type = "Ljung")

    Test de Ljung-Box appliqué à ut pour des lags (de 1 à 30 car au delà on utilise la table de la loi normale)

    ut_TestStat = matrix(0,30) #On génère une matrice pour y stocker les différentes valeurs de Q(m)
for (i in 1:30){
  ut_BoxTest= Box.test(ut, lag = i, type = "Ljung")
  ut_TestStat[i][1] = ut_BoxTest$statistic 
}
plot(ut_TestStat, type = "l", main = "Ljung-Box Test : Ut", xlab = "Lag", ylab = "Statistic X")

    Test de Ljung-Box appliquéà yt pour des lags (de 1 à 30 car au delà on utilise la table de la loi normale)

    yt_TestStat = matrix(0,30) #On génère une matrice pour y stocker les différentes valeurs de Q(m)
for (i in 1:30){
  yt_BoxTest= Box.test(yt, lag = i, type = "Ljung")
  yt_TestStat[i][1] = yt_BoxTest$statistic 
}
plot(yt_TestStat, type = "l", main = "Ljung-Box Test : yt", xlab = "Lag", ylab = "Statistic X")

    • Exercice 2 : Import, tracé et analyse de données.

    Questions :
    1. Import the database contained in the cvs file using the read.csv function. This database is a collection of stock prices.
    2. Plot the third and fifth columns of this dataframe using either the ts_plot or the ggplot2 framework.
    3. Plot the density of the selected stock daily returns ? Using a quantile plot, check the normality of the empirical distribution of the daily returns. Display the four charts within the same output window 6
    4. Calculate the autocovariance and the autocorrelation functions of the daily returns and calculate the Ljung-Box text for the various numbers of lags.
    5. Plot the same ACF of the daily squared returns and calculate the Ljung-Box text for the various numbers of lags.
    6. What are your conclusions regarding the stationary of the daily and squared daily returns

    1 - Import des donnees

    data <- read.csv("data.csv", header = TRUE, sep = ",")
#head(data,5)

    2 - Représentation graphique du cours des actions 3 et 5 ( LVMH & ROG )

    library(ggplot2)
library (dplyr)
library (hrbrthemes)
library(reshape2)
LVMH = data$F.LVMH
ROG = data$S.ROG
dates = as.Date(data$X, format = "%m/%d/%Y")

ggplot(data, aes(x= dates)) + geom_line(aes(y = LVMH, col = "LVMH")) + geom_line(aes(y = ROG, col = "ROG")) + ylab("Stock Price") + ggtitle("StockPrice evolution : LVMH & ROG") + scale_color_manual(values = c('LVMH' = 'darkblue', 'ROG' = 'red')) + labs(color = 'Stock')

    3 - Calcul des rendements moyens et représentations graphiques

    LVMH.Daily_returns <- matrix(0,length(LVMH))
ROG.Daily_returns <- matrix(0,length(ROG))

for (i in 2:length(LVMH)){
  LVMH.Daily_returns[i] = (LVMH[i]-LVMH[i-1])/LVMH[i-1]
}

for (i in 2:length(ROG)){
  ROG.Daily_returns[i] = (ROG[i]-ROG[i-1])/ROG[i-1]
}
 # On determine les densités des rendements : 
p1 <- density(LVMH.Daily_returns)
p2 <- density(ROG.Daily_returns)
    qqnorm() : On appelle QQ-Plot normal le diagramme qui permet de comparer la distribution des données d’un lot à la distribution dite normale ou gaussienne.
  • Dans le cas ci-desosus on compare les deux densités de rendements des actions LVMH & ROG à une densité de loi normale.
    • On obtient alors la représenation graphique suivante :
    par(mfrow = c(2,2))
plot(p1, main = "LVMH returns density", col = "blue")
plot(p2, main = "ROG returns density", col = "green")
qqnorm(LVMH.Daily_returns, ylab = "LVMH Returns quantiles", main = "Normal Q-Q Plot (LVMH)", col = "blue")
qqnorm(ROG.Daily_returns,  ylab = "ROG Returns quantiles", main = "Normal Q-Q Plot (ROG)", col = "green")

    On remarque dans un premier temps que les densités des deux actions sont celles de lois normales. Cette conjecture est confirmée par les QQ-plot.

    4- Autocovariance et Ljung-Box pour différents lags :

    On commence par représenter les ACF respectives des rendements journaliers des actions LVMH et ROG :

    par(mfrow = c(1,2))
acf(LVMH.Daily_returns, lag = 5, main = "LVMH daily returns")
acf(ROG.Daily_returns, lag = 5, main = "ROG daily returns")
    On conjecture que les autocorrélations sont à première vue toutes significativement nulles dès le premier lag.

    Ensuite on applique le test de Ljung-Box pour des lags variants (1~30):

    LB.LVMH_returns = matrix(0,30)
LB.ROG_returns = matrix(0,30)

for (i in 1:30){
  #Evaluation du test de LB pour le lag i.
  LB.boxLvmh = Box.test(LVMH.Daily_returns, lag = i, type = "Ljung")
  LB.boxROG = Box.test (ROG.Daily_returns, lag = i, type = "Ljung")
  
  #Attribution des quantiles : 
  LB.LVMH_returns[i] = LB.boxLvmh$statistic
  LB.ROG_returns[i] = LB.boxROG$statistic
}

    Enfin, on obtient les représentations graphiques suivantes :

    df <- data.frame( x = 1:30, y = LB.LVMH_returns, y2 = LB.ROG_returns)
ggplot(df, aes(x = df$x)) + geom_line(aes(y = df$y,col = "LVMH_statistics")) + geom_line(aes(y= df$y2, col = "ROG_statistics")) + xlab("Lag") + ylab("Statistic X") + ggtitle("Ljung-Box Test ") + scale_color_manual(values = c('LVMH_statistics' = 'darkgreen', 'ROG_statistics' = 'blue')) + labs(color = "Statistics")

    5- Etudes des rendments quadratiques journaliers :

    LVMH.dr.squared = LVMH.Daily_returns^2
acf(LVMH.dr.squared, lag = 10, main = "ACF Quadratic Returns (LVMH)")

    ROG.dr.squared = ROG.Daily_returns^2
acf(ROG.dr.squared, lag = 10, main = "ACF Quadratic Returns (ROG)")

    On peut alors représenter la variation des statistiques relativement à la varaition du lag pour un test de Ljung-Box :

    LB.LVMH_squared_returns = matrix(0,30)
LB.ROG_squared_returns = matrix(0,30)

for (i in 1:30){
  #Evaluation du test de LB pour le lag i.
  LB.boxLvmh = Box.test(LVMH.dr.squared, lag = i, type = "Ljung")
  LB.boxROG = Box.test (ROG.dr.squared, lag = i, type = "Ljung")
  
  #Attribution des quantiles : 
  LB.LVMH_squared_returns[i] = LB.boxLvmh$statistic
  LB.ROG_squared_returns[i] = LB.boxROG$statistic
}

df <- data.frame( x = 1:30, y = LB.LVMH_squared_returns, y2 = LB.ROG_squared_returns)

ggplot(df, aes(x = df$x)) + geom_line(aes(y = df$y,col = "LVMH_statistics")) + geom_line(aes(y= df$y2, col = "ROG_statistics")) + xlab("Lag") + ylab("Statistic X") + ggtitle("Ljung-Box Test \n(squared returns) ") + scale_color_manual(values = c('LVMH_statistics' = 'darkgreen', 'ROG_statistics' = 'blue')) + labs(color = "Statistics")

    • Exercice 3 : EMV

    1- The likelihood of a sample generated by a exponential distribution

    Dans un premier temps, on définit la vraissemblance pour un sample généré par une loi exponentielle. On suppose les variables aléatoires indépendantes et identiquement distribuées.

    #On implémente la fonction à densité exponentielle

densite_expo <- function(x,theta){
  result = theta*exp(-theta*x)
  result
}

# On implémente la fonction de vraissemblance
likelyhood <- function(sample,theta){
  returned = 1
  for (x in sample){
    returned = returned*densite_expo(x,theta)
  }
  returned
  }

    2- the parameter θ :

    le calcul de la vraissemblance pour un echantillon de variables aléatoires indépendantes aboutit à une logvraissemblance :“nlog(θ) - θSomme(xi)” . Donc la dérivée de la logvraissemblance est : “n/θ - Somme(xi). Or on cherche la valeur de théta la plus vraissemblance donc cette dernière corréspond à un extremum de la vraissemblance (ou log(vraissemblance)). On obtient ainsi par simple résolution que l’estimation de” 
    # Or, théta est le paramètre inconnu que l'on doit estimer, on procède alors de la sorte :
  
Estimation_expo <- function(sample){
  xn = sum(sample)/length(sample) #estimateur de la moyenne empirique
  theta = 1/xn
  theta
}

    2- Answer to the same questions considering this time a sample generated by a Gaussian distribution.

    Density/likelyhood :

    #densité : 
densite_gauss <- function(x, nu, sigma){
  densite = (sqrt(2*pi)*sigma)^-1 * exp((1/2)* ((x-nu)/sigma)^2)
}

#Vraissemblance :

likelyhood_gauss <- function(sample, nu, sigma){
  result = 1
  for (x in sample){
    result = result * densite_gauss(x,nu,sigma)
  }
  result
}
    Par la suite, la recherche d’extremums pour les deux paramètres aboutit à :
  • nu : somme(xi)/n
  • sigma : somme((xi-nu)^2)/n

    • On obtient alors :

    nu_sigma <- function(sample){
  if (length(sample) != 0){
    xn = sum(sample)/length(sample)
    sig = 0
    for (x in sample){
      sig = sig + ((x-xn)^2)/length(sample)
    }
  }
  nu = xn
  sigma = sig
  list(nu,sigma)
}